{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## _*Tapering an Operator*_\n",
    "\n",
    "This notebook demonstrates how symmetries can be taken advantage of to reduce the size (number of qubits needed) for an Operator when using Qiskit Chemistry.\n",
    "\n",
    "This notebook has been written to use the PYSCF chemistry driver."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Couldn't find cython int routine\n",
      "Couldn't find cython int routine\n"
     ]
    }
   ],
   "source": [
    "# import common packages\n",
    "import itertools\n",
    "import logging\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "from qiskit import BasicAer\n",
    "\n",
    "from qiskit.aqua import QuantumInstance\n",
    "from qiskit.aqua.operators import Z2Symmetries, WeightedPauliOperator\n",
    "from qiskit.aqua.algorithms import VQE, NumPyMinimumEigensolver\n",
    "from qiskit.aqua.components.optimizers import COBYLA\n",
    "\n",
    "from qiskit.chemistry.drivers import PySCFDriver, UnitsType\n",
    "from qiskit.chemistry.core import Hamiltonian, TransformationType, QubitMappingType \n",
    "from qiskit.chemistry.components.variational_forms import UCCSD\n",
    "from qiskit.chemistry.components.initial_states import HartreeFock"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# using driver to get fermionic Hamiltonian\n",
    "driver = PySCFDriver(atom='Li .0 .0 .0; H .0 .0 1.6', unit=UnitsType.ANGSTROM,\n",
    "                     charge=0, spin=0, basis='sto3g')\n",
    "molecule = driver.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Originally requires 8 qubits\n",
      "Electronic Hamiltonian: Representation: paulis, qubits: 8, size: 276\n"
     ]
    }
   ],
   "source": [
    "core = Hamiltonian(transformation=TransformationType.FULL, qubit_mapping=QubitMappingType.PARITY, \n",
    "                   two_qubit_reduction=True, freeze_core=True)\n",
    "qubit_op, _ = core.run(molecule)\n",
    "\n",
    "print(\"Originally requires {} qubits\".format(qubit_op.num_qubits))\n",
    "print(qubit_op)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find the symmetries of the qubit operator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Z2 symmetries found:\n",
      "ZIZIZIZI\n",
      "ZZIIZZII\n",
      "single qubit operators found:\n",
      "IIIIIIXI\n",
      "IIIIIXII\n",
      "cliffords found:\n",
      "ZIZIZIZI\t(0.7071067811865475+0j)\n",
      "IIIIIIXI\t(0.7071067811865475+0j)\n",
      "\n",
      "ZZIIZZII\t(0.7071067811865475+0j)\n",
      "IIIIIXII\t(0.7071067811865475+0j)\n",
      "\n",
      "single-qubit list: [1, 2]\n"
     ]
    }
   ],
   "source": [
    "z2_symmetries = Z2Symmetries.find_Z2_symmetries(qubit_op)\n",
    "print('Z2 symmetries found:')\n",
    "for symm in z2_symmetries.symmetries:\n",
    "    print(symm.to_label())\n",
    "print('single qubit operators found:')\n",
    "for sq in z2_symmetries.sq_paulis:\n",
    "    print(sq.to_label())\n",
    "print('cliffords found:')\n",
    "for clifford in z2_symmetries.cliffords:\n",
    "    print(clifford.print_details())\n",
    "print('single-qubit list: {}'.format(z2_symmetries.sq_list))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the found symmetries, single qubit operators, and cliffords to taper qubits from the original qubit operator. For each Z2 symmetry one can taper one qubit. However, different tapered operators can be built, corresponding to different symmetry sectors. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of qubits of tapered qubit operator: 6\n",
      "Number of qubits of tapered qubit operator: 6\n",
      "Number of qubits of tapered qubit operator: 6\n",
      "Number of qubits of tapered qubit operator: 6\n"
     ]
    }
   ],
   "source": [
    "tapered_ops = z2_symmetries.taper(qubit_op)\n",
    "for tapered_op in tapered_ops:\n",
    "    print(\"Number of qubits of tapered qubit operator: {}\".format(tapered_op.num_qubits))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The user has to specify the symmetry sector he is interested in. Since we are interested in finding the ground state here, let us get the original ground state energy as a reference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=== GROUND STATE ENERGY ===\n",
      " \n",
      "* Electronic ground state energy (Hartree): -8.874303870396\n",
      "  - computed part:      -1.078084301625\n",
      "  - frozen energy part: -7.796219568771\n",
      "  - particle hole part: 0.0\n",
      "~ Nuclear repulsion energy (Hartree): 0.992207270475\n",
      "> Total ground state energy (Hartree): -7.882096599921\n"
     ]
    }
   ],
   "source": [
    "ee = NumPyMinimumEigensolver(qubit_op)\n",
    "result = core.process_algorithm_result(ee.run())\n",
    "print(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let us iterate through all tapered qubit operators to find out the one whose ground state energy matches the original (un-tapered) one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Lowest eigenvalue of the 0-th tapered operator (computed part) is -1.078084301625\n",
      "Lowest eigenvalue of the 1-th tapered operator (computed part) is -0.509523578167\n",
      "Lowest eigenvalue of the 2-th tapered operator (computed part) is -0.912078232998\n",
      "Lowest eigenvalue of the 3-th tapered operator (computed part) is -0.912078232998\n",
      "The 0-th tapered operator matches original ground state energy, with corresponding symmetry sector of [1, 1]\n"
     ]
    }
   ],
   "source": [
    "smallest_eig_value = 99999999999999\n",
    "smallest_idx = -1\n",
    "for idx in range(len(tapered_ops)):\n",
    "    ee = NumPyMinimumEigensolver(tapered_ops[idx])\n",
    "    curr_value = ee.run().eigenvalue.real\n",
    "    if curr_value < smallest_eig_value:\n",
    "        smallest_eig_value = curr_value\n",
    "        smallest_idx = idx\n",
    "    print(\"Lowest eigenvalue of the {}-th tapered operator (computed part) is {:.12f}\".format(idx, curr_value))\n",
    "    \n",
    "the_tapered_op = tapered_ops[smallest_idx]\n",
    "the_coeff = tapered_ops[smallest_idx].z2_symmetries.tapering_values\n",
    "print(\"The {}-th tapered operator matches original ground state energy, with corresponding symmetry sector of {}\".format(smallest_idx, the_coeff))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, one can run multiple VQE instances to find the lowest eigenvalue sector. \n",
    "Here we just validate that `the_tapered_op` reach the smallest eigenvalue in one VQE execution with the UCCSD variational form, modified to take into account of the tapered symmetries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# setup initial state\n",
    "init_state = HartreeFock(num_orbitals=core._molecule_info['num_orbitals'],\n",
    "                    qubit_mapping=core._qubit_mapping, two_qubit_reduction=core._two_qubit_reduction,\n",
    "                    num_particles=core._molecule_info['num_particles'], sq_list=the_tapered_op.z2_symmetries.sq_list)\n",
    "\n",
    "# setup variationl form\n",
    "var_form = UCCSD(num_orbitals=core._molecule_info['num_orbitals'], \n",
    "                 num_particles=core._molecule_info['num_particles'],\n",
    "                 active_occupied=None, active_unoccupied=None, initial_state=init_state,\n",
    "                 qubit_mapping=core._qubit_mapping, two_qubit_reduction=core._two_qubit_reduction, \n",
    "                 num_time_slices=1, z2_symmetries=the_tapered_op.z2_symmetries)\n",
    "\n",
    "# setup optimizer\n",
    "optimizer = COBYLA(maxiter=1000)\n",
    "\n",
    "# set vqe\n",
    "algo = VQE(the_tapered_op, var_form, optimizer)\n",
    "\n",
    "# setup backend\n",
    "backend = BasicAer.get_backend('statevector_simulator')\n",
    "quantum_instance = QuantumInstance(backend=backend)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "algo_result = algo.run(quantum_instance)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=== GROUND STATE ENERGY ===\n",
      " \n",
      "* Electronic ground state energy (Hartree): -8.874303841496\n",
      "  - computed part:      -1.078084272725\n",
      "  - frozen energy part: -7.796219568771\n",
      "  - particle hole part: 0.0\n",
      "~ Nuclear repulsion energy (Hartree): 0.992207270475\n",
      "> Total ground state energy (Hartree): -7.882096571021\n",
      "The parameters for UCCSD are:\n",
      "[ 0.03803094  0.00360661  0.0382837   0.00369849 -0.03608325  0.05942172\n",
      " -0.02729715 -0.02732059  0.05965191 -0.11498381]\n"
     ]
    }
   ],
   "source": [
    "result = core.process_algorithm_result(algo_result)\n",
    "print(result)\n",
    "\n",
    "print(\"The parameters for UCCSD are:\\n{}\".format(algo_result.optimal_point))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
